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We  present  a  platform  for  the  simulation  of  quantum  magnetism 
with  full  control  of  interactions  between  pairs  of  spins  at  arbitrary 
distances  in  ID  and  2D  lattices.  In  our  scheme,  two  internal  atomic 
states  represent  a  pseudospin  for  atoms  trapped  within  a  photonic 
crystal  waveguide  (PCW).  With  the  atomic  transition  frequency 
aligned  inside  a  band  gap  of  the  PCW,  virtual  photons  mediate 
coherent  spin-spin  interactions  between  lattice  sites.  To  obtain  full 
control  of  interaction  coefficients  at  arbitrary  atom-atom  separa¬ 
tions,  ground-state  energy  shifts  are  introduced  as  a  function  of 
distance  across  the  PCW.  In  conjunction  with  auxiliary  pump  fields, 
spin-exchange  versus  atom-atom  separation  can  be  engineered 
with  arbitrary  magnitude  and  phase,  and  arranged  to  introduce 
nontrivial  Berry  phases  in  the  spin  lattice,  thus  opening  new  ave¬ 
nues  for  realizing  topological  spin  models.  We  illustrate  the  broad 
applicability  of  our  scheme  by  explicit  construction  for  several  well- 
known  spin  models. 

nanophotonics  |  quantum  matter  |  cold  atoms  |  quantum  many-body  | 
quantum  spin 

Quantum  simulation  has  become  an  important  theme  for  re¬ 
search  in  contemporary  physics  (1).  A  quantum  simulator 
consists  of  quantum  particles  (e.g.,  neutral  atoms)  that  interact  by 
way  of  a  variety  of  processes,  such  as  atomic  collisions.  Such  pro¬ 
cesses  typically  lead  to  short-range,  nearest-neighbor  interactions 
(2-6).  Alternative  approaches  for  quantum  simulation  use  dipolar 
quantum  gases  (7,  8),  polar  molecules  (9-11),  and  Rydberg  atoms 
(12-15),  leading  to  interactions  that  typically  scale  as  1/r3,  where  r  is 
the  interparticle  separation.  For  trapped  ion  quantum  simulators 
(16-20),  tunability  in  a  power  law  scaling  of  r~t]  with  0<rj<3  can  in 
principle  be  achieved.  Beyond  simple  power  law  scaling,  it  is  also 
possible  to  engineer  arbitrary  long-range  interactions  mediated  by 
the  collective  phonon  modes,  which  can  be  achieved  by  independent 
Raman  addressing  on  individual  ions  (21). 

Using  photons  to  mediate  controllable  long-range  interactions 
between  isolated  quantum  systems  presents  yet  another  approach 
for  assembling  quantum  simulators  (22).  Recent  successful  ap¬ 
proaches  include  coupling  ultracold  atoms  to  a  driven  photonic 
mode  in  a  conventional  mirror  cavity,  thereby  creating  quantum 
many-body  models  (using  atomic  external  degrees  of  freedom) 
with  cavity-field-mediated  infinite-range  interactions  (23). 
Finite-range  and  spatially  disordered  interactions  can  be  realized 
by  using  multimode  cavities  (24).  Recent  demonstrations  on 
coupling  cold  atoms  to  guided  mode  photons  in  photonic  crystal 
waveguides  (25,  26)  and  cavities  (27,  28)  present  promising  ave¬ 
nues  (using  atomic  internal  degrees  of  freedom)  due  to  unprec¬ 
edented  strong  single  atom-photon  coupling  rate  and  scalability. 
Related  efforts  also  exists  for  coupling  solid-state  quantum  emit¬ 
ters,  such  as  quantum  dots  (29,  30)  and  diamond  nitrogen-vacancy 
centers  (31,  32),  to  photonic  crystals.  Scaling  to  a  many-body 
quantum  simulator  based  on  solid-state  systems,  however,  still 
remains  elusive.  Successful  implementations  can  be  found  in 
the  microwave  domain,  where  superconducting  qubits  behave  as 
artificial  atoms  strongly  coupled  to  microwave  photons  propa¬ 
gating  in  a  network  formed  by  superconducting  resonators  and 
transmission  lines  (33-35). 


Here,  we  propose  and  analyze  a  physical  platform  for  simulating 
long-range  quantum  magnetism  in  which  full  control  is  achieved  for 
the  spin-exchange  coefficient  between  a  pair  of  spins  at  arbitrary 
distances  in  ID  and  2D  lattices.  The  enabling  platform,  as  described 
in  refs.  36  and  37,  is  trapped  atoms  within  photonic  crystal  wave¬ 
guides  (PCWs),  with  atom-atom  interactions  mediated  by  photons 
of  the  guided  modes  (GMs)  in  the  PCWs.  As  illustrated  in  Fig.  1 A 
and  B,  single  atoms  are  localized  within  unit  cells  of  the  PCWs  in  ID 
and  2D  periodic  dielectric  structures.  At  each  site,  two  internal 
atomic  states  are  treated  as  pseudospin  states,  with  spin- 1/2  con¬ 
sidered  here  for  definiteness  (e.g.,  states  | g)  and  | s)  in  Fig.  1C). 

Our  scheme  uses  strong,  and  coherent  atom-photon  interactions 
inside  a  photonic  band  gap  (36-40),  and  long-range  transport 
property  of  GM  photons  for  the  exploration  of  a  large  class  of 
quantum  magnetism.  This  is  contrary  to  conventional  hybrid 
schemes  based  on,  for  example,  arrays  of  high  finesse  cavities 
(41-44)  in  which  the  pseudospin  acquires  only  the  nearest  (or  at 
most  the  next-nearest)  neighbor  interactions  due  to  strong  expo¬ 
nential  suppression  of  photonic  wave  packet  beyond  single  cavities. 

In  its  original  form  (36-40),  the  localization  of  pseudospin  is 
effectively  controlled  by  single-atom  defect  cavities  (36).  The  cavity 
mode  function  can  be  adjusted  to  extend  over  long  distances  within 
the  PCWs,  thereby  permitting  long-range  spin  exchange  interac¬ 
tions.  The  interaction  can  also  be  tuned  dynamically,  via  external 
addressing  beams,  to  induce  complex  long-range  spin  transport, 
which  we  describe  in  the  following  (36,  37). 

To  engineer  tunable,  long-range  spin  Hamiltonians,  we  use  an 
atomic  A  scheme  and  two-photon  Raman  transitions,  where  an  atom 
flips  its  spin  state  by  scattering  one  photon  from  an  external  pump 
field  into  the  GMs  of  a  PCW.  The  GM  photon  then  propagates 
within  the  waveguide,  inducing  spin  flip  in  an  atom  located  at  a 


Significance 

Cold  atoms  trapped  along  a  photonic  crystal  waveguide  can  be 
used  to  simulate  long-range  quantum  magnetism  with  pair¬ 
wise-tunable  spin-spin  interactions  mediated  by  guided  virtual 
photons  in  a  photonic  band  gap.  Using  a  two-photon  Raman 
addressing  scheme,  the  proposed  atom-nanophotonic  system 
can  achieve  arbitrary  and  dynamic  control  on  the  strength, 
phase,  and  length  scale  of  spin  interactions.  This  promises  new 
avenues  for  engineering  a  large  class  of  spin  Hamiltonians, 
including  those  exhibiting  topological  order  or  frustrated  long- 
range  magnetism. 
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Fig.  1.  Photon-mediated  atom-atom  interactions  in  (A)  ID  and  (B)  2D  PGA/s. 
(C)  Atomic-level  scheme:  atomic  dipole  |s)  <-»•  |e)  is  coupled  to  an  external  pump, 
| g)  ^  |e)  coupled  to  a  GM  photon,  and  r*,  the  excited  state  decay  rate  to  free  space 
and  leaky  modest  (D)  Simplified  band  structure  rw  (k)  near  the  band  edge  k  =  kc  and 
&>(kc)  =a>c-  Atomic  transition  frequency  (oeg  =  (oe-(og  lies  within  the  band  gap. 


distant  site  via  the  reverse  two-photon  Raman  process.  When  we 
align  the  atomic  resonant  frequency  inside  the  photonic  band  gap,  as 
depicted  in  Fig.  ID,  only  virtual  photons  can  mediate  this  remote 
spin  exchange  and  the  GM  dynamics  are  fully  coherent,  effectively 
creating  a  spin  Hamiltonian  with  long-range  interactions.  As  dis¬ 
cussed  in  refs.  36  and  37,  the  overall  strength  and  length  scale  of  the 
spin-exchange  coefficients  can  be  tuned  by  an  external  pump  field, 
albeit  within  the  constraints  set  by  a  functional  form  that  depends  on 
the  dimensionality  and  the  photonic  band  structure.  These  con¬ 
straints  may  limit  our  ability  to  explore  novel  quantum  phases  and 
nonequilibrium  dynamics  in  various  spin  models,  because  many  ef¬ 
fects  display  strong  dependencies  on  the  functional  form  of  long- 
range  interactions  (45-50).  It  is  therefore  highly  desirable  to  obtain 
full  control  of  interactions  without  the  need  to  investigate  over  a 
wide  range  of  PCW  designs  with  different  photonic  band  structures. 

To  fully  control  spin-exchange  coefficients  at  arbitrary  separations, 
here  we  adopt  a  Raman-addressing  scheme  similarly  discussed  for 
cold  atoms  and  trapped  ions  (51-55).  We  introduce  atomic  ground- 
state  energy  shifts  as  a  function  of  distance  across  the  PCW.  Due  to 
conservation  of  energy,  these  shifts  suppress  reverse  two-photon 
Raman  processes  in  the  original  scheme  (36,  37),  forbidding  spin 
exchange  within  the  entire  PCW.  However,  we  can  selectively  acti¬ 
vate  certain  spin-exchange  interactions  J(rm >n)  between  atom  pairs 
(m,n)  separated  by  rm,n,  by  applying  an  auxiliary  sideband  whose 
frequency  matches  that  of  the  original  pump  plus  the  ground-state 
energy  shift  between  the  atom  pairs.  This  allows  us  to  build  a  pre¬ 
scribed  spin  Hamiltonian  with  interaction  terms  “one  by  one.”  Note 
that  each  sideband  in  a  Raman-addressing  beam  can  be  easily  in¬ 
troduced,  for  example,  by  an  electro-optical  modulator.  By  in¬ 
troducing  multiple  sidebands  and  by  controlling  their  frequencies, 
amplitudes,  and  relative  phases,  we  can  engineer  spin  Hamiltonians 
with  arbitrary,  complex  interaction  coefficients  J(rm^).  Depending 
on  the  dimensionality  and  the  type  of  spin  Hamiltonians,  our  scheme 
requires  only  one  or  a  few  Raman  beams  to  generate  the  desired 
interactions.  Furthermore,  by  properly  choosing  the  propagation 
phases  of  the  Raman  beams,  we  can  imprint  geometric  phases  in  the 
spin  system,  thus  providing  unique  opportunities  for  realizing  topo¬ 
logical  spin  models. 

We  substantiate  the  broad  applicability  of  our  methods  by 
explicit  elaboration  of  the  set  of  pump  fields  required  to  realize 
well-known  spin  Hamiltonians.  For  ID  spin  chains,  we  consider 
the  implementation  of  the  Haldane-Shastry  model  (56,  57).  For 
2D  spin  lattices,  we  elaborate  the  configurations  for  realizing 
topological  flat  bands  (58,  59)  in  Haldane’s  spin  model  (56),  as 
well  as  a  “checkerboard”  chiral-flux  lattice  (58,  59).  We  also 
consider  a  2D  XXZ  spin  Hamiltonian  with  J(rm,n)  <x  l/rjnn  and 


rj  =  1,2,3  (60).  In  addition,  we  report  numerical  results  on  the  rj 
dependence  of  its  magnetization  diagram. 

Controlling  Spin-Spin  Interaction  Through  Multifrequency 
Driving 

In  the  following,  we  discuss  how  to  achieve  full  control  of  interac¬ 
tions  by  multifrequency  pump  fields.  We  assume  (i)  N  atoms  trapped 
in  either  a  ID  or  2D  PCW,  as  depicted  in  Fig.  1  A  and  B ,  with  a 
spatially  dependent  ground-state  energy  shift  cog.  For  simplicity,  we 
assume  one  atom  per  unit  cell  of  the  PCW,  although  this  assumption 
can  be  relaxed  afterward;  (ii)  the  structure  is  engineered  (22-28) 
such  that  the  GM  polarization  is  coupled  to  the  atomic  dipole, 

| g)  ++  \e),  as  shown  in  Fig.  1C,  and,  under  rotating  wave  approxi¬ 
mation,  is  described  by  the  following  Hamiltonian  (using  h  =  1): 

Him  =  y]gk(rn)«k<4  +  h.c.,  [1] 

k  ,n 

where  gk(rn)  =gkelk  Tn  is  the  single-photon  coupling  constant  at 
site  location  rn,  with  n  being  the  site  index;  ak,  the  GM  field 
operator;  and  o^b  =  \a)n(b\,  the  atomic  operators  with  a,b  being 
one  of  the  g,s,e  states.  Moreover,  as  in  refs.  36  and  37,  we 
assume  (ui)  there  is  another  hyperfine  level  | s),  addressed  by  a 
Raman  field  with  coupling  strength  Q  as  follows: 

^dW=E(£y^/w+hx)  [2] 

where  col  is  the  main  driving  frequency.  The  Raman  field  Q(t) 
contains  mp  frequency  components  that  are  introduced  to 
achieve  full  control  of  the  final  effective  spin  Hamiltonian.  Full 
dependence  of  Q(f)  can  be  written  as  follows: 

Hlp~  1 

0(0  =  E  [3] 

a— 0 

where  cba  are  the  detunings  of  the  sidebands  from  the  main  frequency 
coL  such  that  coq  =  0,  and  Qa,  the  complex  amplitudes. 

We  can  adiabatically  eliminate  the  excited  states  \e)  and  the 
photonic  GMs  under  the  condition  that  (iv)  max{|Q|,  \coa  -  d)p\}  <C 
|  A |  =  \a>e  -  &>l|-  This  condition  guarantees  that,  first,  the  excited 
state  is  only  virtually  populated,  and  that,  second,  the  time  de¬ 
pendence  induced  by  the  sideband  driving  is  approximately  constant 
over  the  timescale  A-1.  As  discussed  in  refs.  36  and  37,  if  cop  -  mg 
lies  in  the  photonic  band  gap,  photon-mediated  interactions  by  GMs 
are  purely  coherent^  Under  the  Bom-Markov  approximation,  we 
then  arrive  at  an  effective  AT  Hamiltonian  (SI  Appendix  A:  Complete 
Derivation  of  Final  Time-Dependent  Hamiltonian): 

N  mP- 1 

Hxr(t)=  [4] 

m,  n^m  a,  /i=0 

where  we  have  defined  Xa  =  Qa/(2A);  cog^n=oog( rn)  is  the  site- 
dependent  ground-state  energy  shift,  and  J(rm/l)  is  the  atom- 
GM  photon  coupling  strength  (36,  37)  that  typically  depends 
on  atomic  separation  rm  „  =rm-rn. 

We  focus  on  “sideband  engineering”  and  treat  J(rm/l)  as  ap¬ 
proximately  constant  over  atomic  separations  considered.* *  This 
is  valid  as  long  as  the  farthest  atomic  separation  with  nonzero 
engineered  interaction  is  much  smaller  than  the  decay  length 


+To  simplify  the  discussion,  in  this  paper,  we  neglect  decoherence  effects  caused  by 
atomic  emission  into  free  space  and  leaky  modes  as  well  as  photon  loss  due  to  imper¬ 
fections  in  the  PCW.  These  effects  were  both  carefully  discussed  in  refs.  36  and  37, 
suggesting  the  number  of  spin-exchange  cycles  in  the  presence  of  decoherence  can  re¬ 
alistically  reach  yf’»35~100  using  ultra-high  Q  PCWs. 

*One  may  also  replace  a  PCW  with  a  single-mode  nanophotonic  cavity,  operating  in  the 
strong  dispersive  regime  (61,  62),  to  achieve  constant  GM  coupling  J  independent  of 
|rm#n|.  Realistic  nanophotonic  cavity  implementations  will  be  considered  elsewhere. 
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[8] 


\e) 


resonant  Raman-scattering  processes.  Spin  exchanges  (/A)  \sn,gm)  -*■  \9n,sm) 
and  ( B )  \gn,sm)  |sn,gm)  are  allowed  only  when  the  condition  -<ygn  = 
cbp-cba  is  satisfied.  Qa/ A  and  Q.p/ A  control  the  exchange  rate. 


m.p-1 

Jm,n  —  ^  ^  XaX pi 8 n  —m ,fi—a  • 
a,f}= 0 

It  can  be  shown  that  the  set  of  equations  defined  by  Eq.  8  has  at  least 
one  solution  for  any  arbitrary  choice  of  Jms,  that  is,  by  choosing 
£20  >  o  and  Jm,n  «  (XoX*_m  +XgXm_n)/.  More  solutions  can 
be  found  by  directly  solving  the  set  of  nonlinear  equations  Eq.  8. 

It  is  important  to  highlight  that  multifrequency  driving  also 
enables  the  possibility  to  engineer  geometrical  phases  and, 
therefore,  topological  spin  models.  If  the  pump  field  propagation 
is  not  perfectly  transverse,  that  is,  kL  •  rm(n)  ^  0  (kL  being  the 
wave  vector  of  the  Raman  field),  the  effective  Hamiltonian  Eq.  7 
acquires  spatial-dependent,  complex  spin-exchange  coefficients 
via  the  phase  of  XaX ^  in  Eq.  8;  see  later  discussions. 

Beyond  an  ideal  setting,  we  now  stress  a  few  potential  error 
sources.  First,  for  practical  situations,  the  gradient  per  site  8  will 
be  a  limited  resource,  making  Eq.  7  not  an  ideal  approximation. 
Careful  Floquet  analysis  on  time-dependent  Hamiltonian  in  Eqs. 
5  and  6  is  required,  to  be  discussed  later.  Second,  there  is  an 
additional  Stark  shift  on  state  | s)  due  to  the  Raman  fields: 


scale  f  of  the  coupling  strength  Here,  A  is  the 

band  curvature  (Fig.  ID),  Ac  =  max{&>c  —  (cop  —  <%«)}  is  the 
maximal  detuning  of  the  band  edge  to  the  frequency  of  coupled 
virtual  photons  that  mediate  interactions  (Fig.  1C),  and  we  have 
assumed  that  the  variation  of  ground-state  energies  cog,n  are  small 
compared  with  Ac.  Exact  functional  form  of  J(rm,n)  can  be  found 
in  refs.  36  and  37,  and  in  SI  Appendix  A:  Complete  Derivation  of 
Final  Time-Dependent  Hamiltonian. 

The  time  dependence  in  Eq.  4  can  be  further  engineered  and 
simplified.  We  note  that  the  interaction  between  two  atoms  n  and  m 
will  be  highly  dependent  on  the  resonant  condition  oog,m  -  cog,n  = 
(bp  -  ooa,  provided  the  ground-state  energy  difference  \cog,n  -  cog/n  |  is 
much  larger  than  the  characteristic  timescale  of  interactions 
\XaX*J\.  The  intuitive  picture  is  depicted  in  Fig.  24:  the  atom 
n  scatters  from  sideband  a  a  photon  with  energy  coL  +  a?a-  oog,n  into 
the  GMs.  When  this  GM  photon  propagates  to  the  atom  m,  it  will 
only  be  rescattered  into  a  sideband  f  that  satisfies  cop  +  8oa  - 
a>g,n  =  cop  +  oop-  cQg,m,  whereas  the  rest  of  the  sidebands  remain 
off-resonant.  Fig.  2 B  depicts  a  reversed  process. 

For  concreteness,  we  discuss  a  ID  case  where  we  assume  (v)  a 
linear  gradient  in  the  ground-state  energy  cog/l  =  n8,  with  8  being 
the  energy  difference  between  adjacent  sites.  The  sidebands  will 
be  chosen  accordingly  such  that  cba  =  a8,  with  aEZ. 

Summing  up,  with  all  these  assumptions  (i-v),  the  resulting 
effective  Hamiltonian  Eq.  4  can  finally  be  rewritten  as  follows: 

Hxv(t)  =  YHxy,p  eipSt>  [5] 

p 

where  Hxy,p  is  the  contribution  that  oscillates  with  frequency  p8. 
Written  explicitly, 


N  rap- 1 

Hxy,p  =  y!  XaXpJ8n -mj-a-p (/^s (Asg .  [6] 

m,n^m  a, (3=0 


In  an  ideal  situation,  the  gradient  per  site  satisfies  8  »  \xax;i\ 
such  that  the  contributions  from  Hxy,p  V/?/0  can  be  neglected. 
Under  these  assumptions,  we  arrive  at  an  effective  time-indepen¬ 
dent  Hamiltonian: 


N 

Hxy{t)  ~Hxy,o  =  ^2  JmsiO&Osg’  [7] 

m,n^m 


where  couplings  Jm ^  can  be  tuned  by  adjusting  the  amplitudes 
and  phases  of  the  sidebands  Xa  as  they  are  given  by  the  following: 


8oos{t)  =  - 


rap-l 

E 


a= 0 


|£4|2 

4A 


rap-1 

E* 


a>(3 


^‘a^‘Pj(wa-Wg)t 

2A 


[9] 


where  9i[.]  indicates  real  part.  We  note  that  the  time-independent 
contribution  in  Eq.  9  can  be  absorbed  into  the  energy  of  cos  without 
significant  contribution  to  the  dynamics,  whereas  the  time-depen¬ 
dent  terms  may  be  averaged  out  over  the  atomic  timescales  that  we 
are  interested  in.  We  will  present  strategies  for  optimizing  the 
choice  of  8,  and  minimizing  detrimental  effects  due  to  undesired 
time-dependent  terms  in  Eqs.  5  and  9  in  later  discussions. 

Independent  Control  of  XX  and  YY  Interactions.  So  far,  we  can  fully 
engineer  an  XY  Hamiltonian  with  equal  weight  between  XX  and 
YY  terms  by  defining  the  Pauli  operators  (<jx,  oy,  of)  =  (<7sg  + 
<7gs,  i(osg  -  <7gs  ),  Ogg  -  <7SS).  We  now  show  flexible  control  of  AX  and 
YY  interactions  with  slight  modifications  in  the  atomic  level  structure 
and  the  Raman-addressing  scheme.  In  particular,  we  use  a  butterfly¬ 
like  level  structure  where  there  are  two  transitions,  \g)  \e)  and 

| s)  \e),  coupled  to  the  same  GM,  as  depicted  in  Fig.  3.  We  will  use 
two  multifrequency  Raman  pump  fields,  Qg(t)  and  Qs  (t),  to  induce 
| g)  ++  \e)  <-»  | s)  and  | s)  \e)  |g)  two-photon  Raman  transitions, 
respectively. 

For  example,  to  control  XX  or  YY  interactions,  we  require  that 
the  two  pump  fields  induce  spin  flips  with  equal  amplitude,  that 
is,  ogs  ±  osg.  This  is  possible  if  we  choose  the  main  frequencies  of 
the  pumps  ( ooL,g  and  copf)  such  that  6l>l,s  =  a>Lg  +  2cog,  and  match 
their  amplitudes  such  that  |^g;Cr|/Ag  =  |^5;a|/A5,  where  A5  = 
coe  —  copj,  Ag  =  cog  —  (copg  +  and  |  A^j  | | • 

Adiabatically  eliminating  the  excited  states  as  well  as  the  GMs, 
we  arrive  at  the  following  Hamiltonian: 

HXX'YY#  =  jf  [jm,n  (og  +  ^4)  (og  +  e-^4)  +  h.c.  ], 

m,n>m 

[10] 

where  is  the  relative  phase  between  the  pumps  fields  Assum¬ 
ing  the  laser  beams  that  generate  the  Raman  fields  are  copropagating 
or  are  both  illuminating  the  atoms  transversely,  that  is,  kL  •  r m,n  =  0, 
we  can  generate  either  X  or  Y  components,  (ofg  ±<t^),  by  setting  the 
phase  (pgs  =  0  or  n;  more  exotic  combinations  are  available  with 
generic  choice  of  (pgs .  Moreover,  if  the  laser  beams  are  not  cop¬ 
ropagating,  they  create  spatially  dependent  phases  (pgs  m .  This  can 
create  site-dependent  XX,  YY,  or  XF  terms. 
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Fig.  3.  Atomic  "butterfly"  level  structure.  Two  pump  fields  Qs  and 
tuned  to  couple  to  the  same  GM  photon,  are  introduced  to  control  XX  and 
YY  interactions  independently. 


Independent  Control  of  ZZ  Interactions.  An  independently  con¬ 
trolled  ZZ  Hamiltonian,  in  combination  with  arbitrary  XY  terms, 
would  allow  us  to  engineer  SU(2)-invariant  spin  models  as  well 
as  a  large  class  of  XXZ  models,  that  is,  the  following: 

HXXZ=HXY+HZZ=  jr  [(2C^<  +  hx)+4X^]- 

m,n>m 

[11] 


In  refs.  36  and  37,  it  was  shown  that  ZZ  interaction  can  be  created 
by  adding  an  extra  pump  field  to  the  | g)  ++  \e)  transition  in  Fig.  1C. 
However,  as  ZZ  terms  in  this  scheme  (36,  37)  do  not  involve  flipping 
atomic  states,  it  is  not  directly  applicable  to  our  multifrequency  pump 
method.  Nonetheless,  because  we  can  generate  XX  and  YY  inter¬ 
actions  independently,  a  straightforward  scheme  to  engineer  Hzz  is 
to  use  single  qubit  rotations  to  rotate  the  spin  coordinates  X++Z  or 
y^Z,  followed  by  stroboscopic  evolutions  (63)  to  engineer  the  full- 
spin  Hamiltonian.  Spin-rotation  can  be  realized,  for  example,  with 
a  collective  microwave  driving  Hmw  =  J]((£2mw/2)e^  +  h.c.),  in 

n 

which  a  ^/2-microwave  pulse  rotates  the  basis  |s)„j 

{(\g)n  +  \s)n)/^Hg)n  +  \s)n)/^}- 
Thus,  an  Hxxz  Hamiltonian  can  be  simulated  using  the  fol¬ 
lowing  stroboscopic  evolution:  {Hxy,  Hzz,  Hxy,  Hzz,  •  •  •  }  in  Nt 
steps  as  schematically  depicted  in  Fig.  4.  As  shown  in  SI  Appendix 
A:  Complete  Derivation  of  Final  Time-Dependent  Hamiltonian , 
the  error  accumulated  in  these  Nt  steps  can  be  bounded  by  the 
following: 


E2< 


N(RJt)2 


N, 


[12] 


where  J  =  max[/m  n]  is  the  largest  energy  scale  of  the  Hamiltonian  we 
want  to  simulate,  andR  is  the  approximate  number  of  atoms  coupled 
through  the  interaction.  For  example,  if  JmA  is  a  nearest-neighbor 
interaction,  R  =  1.  If  Jm,n  oA/\m  -n\\  then  R  cx  Yln= l  V \n which 
typically  grows  much  slower  than  N.  Because  E2ocl  /Nt,  the  Trotter 
error  in  Nt  steps  can  in  principle  be  decreased  to  a  given  accuracy  s 
by  using  enough  steps,  that  is,  Nt  >  ( N(RJt)2)/s . 

More  complicated  stroboscopic  evolutions  may  lead  to  a  more 
favorable  error  scaling  (64-66),  although  in  real  experiments 
there  will  be  a  trade-off  between  minimizing  the  Trotter  error 
and  the  fidelity  of  the  individual  operations  to  achieve  Hxy  and 
Hzz •  As  this  will  depend  on  the  particular  experimental  setup,  we 
will  leave  such  analysis  out  of  current  discussions.  For  illustra¬ 
tion,  we  will  only  consider  the  simplest  kind  of  stroboscopic 
evolution  that  we  depicted  in  Fig.  4. 


Engineering  Spin  Hamiltonians  for  ID  Systems:  The 
Haldane-Shastry  5  =  1/2  Spin  Chain 

In  the  first  example,  we  engineer  a  Haldane-Shastry  spin 
Hamiltonian  in  one  dimension  (56,  57): 


N-lN-m 

whs = z  z j"  [2  ( c?g(%+n  +  h.c.)  +  ,  [13] 

m— 1  n- 1 

where  Jn  =Jo/sm1(nn/N ),  Jq=Jji2  /TV2,  and  N  is  the  number  of 
spins.  The  interaction  strength  decays  slowly  with  approximately 
a  1/r2  dependence  while  satisfying  a  periodic  boundary  condi¬ 
tion.  Such  a  spin  Hamiltonian  is  difficult  to  realize  in  most  phys¬ 
ical  setups  that  interact,  for  example,  via  dipolar  interactions. 

We  can  engineer  the  periodic  boundary  condition  and  the 
long-range  interaction  Jn  directly  using  a  linear  array  of  trapped 
atoms  coupled  to  a  PCW.  To  achieve  this,  we  induce  atomic 
ground-state  energy  shift  m8  according  to  the  spin  index  m ,  and 
then  uniformly  illuminate  the  trapped  atoms  with  an  external 
pump  consisting  of  N  frequency  components  cba  =  ad ,  each  with 
an  amplitude  denoted  by  Qa  and  a  =  0,1,  Regardless 

of  the  position  of  atoms,  all  pump  pairs  with  frequency  dif¬ 
ference  nS  contribute  to  the  spin  interaction  Jn.  Considering 
first  the  XY  terms,  and  according  to  Eq.  7,  we  demand  the 
following: 


N-n-l 


jn*j  ]T  xax*+n  =— 


Jo 


a— 0 


sin2(mr/A/y 


[14] 


where  J  is  the  GM  photon  coupling  rate  (Eq.  8)  that  we  will 
assume  to  be  a  constant  for  the  simplicity  of  discussions.  This 
requires  that  the  physical  size  of  the  spin  chain  be  small 
compared  with  the  decay  length  of  J.  That  is,  Nd  <C  £  where 
d  is  the  atomic  separation.  It  is  then  straightforward  to  find 
the  required  pump  amplitudes  Qa  (or  equivalently  Xa )  by 
solving  Eq.  14  for  all  n.  Notice  that  the  system  of  equations 
Eq.  14  is  overdetermined,  and  therefore  one  can  find  several 
solutions  of  it.  However,  we  choose  the  solution  that  mini¬ 
mizes  the  total  intensity  |^«|2-  Fig-  5  shows  that  the  total  in¬ 
tensity  converges  to  a  constant  value  for  large  N,  as  a  result  of 
decreasing  sideband  amplitudes  for  decreasing  1/r2  interaction 
strengths.  This  is  confirmed  in  Fig.  5  as  we  see  the  growth  of 
the  ratio  between  maximum  and  minimum  sideband  amplitudes 
when  N  increases.  The  same  external  pump  configuration  can  also 
be  used  to  induce  the  ZZ  terms  by  applying  stroboscopic  proce¬ 
dures  as  discussed  in  the  previous  section. 

Engineering  Spin  Hamiltonians  for  2D  Systems:  Topological 
and  Frustrated  Hamiltonians 

In  the  following,  we  discuss  specific  examples  for  engineering  2D 
spin  Hamiltonians  that  are  topologically  nontrivial.  In  particular, 
we  discuss  two  chiral-flux  lattice  models  that  require  long-range 


Nt 

Hxy 

Hzz 

Hxy 

Hzz 

Hxy 

Hzz 

Time  t 

Fig.  4.  Scheme  for  generating  an  XXZ  spin  Hamiltonian  using  a  stroboscopic 
evolution.  The  scheme  contains  periodic  applications  of  a  multifrequency  Raman 
field  to  induce  the  HXy  interaction  (in  green),  two  fast  microwave  pulses  (or 
optical  two-photon  transition)  forming  Hmw  that  uniformly  rotate  the  spin  basis 
{|9>„.  |s)„}  « {(|g)„  +  |s)„)/n/2,  (Ig)„  -  |s)„)/v5}  back  and  forth  (in  blue),  and  a 
butterfly-like  pumping  scheme  that  applies  HZz  in  the  rotated  basis. 
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hopping  terms  to  engineer  single  particle  flat-bands  with  nonzero 
Chern  numbers,  which  are  key  ingredients  to  realizing  fractional 
quantum  Hall  effects  (FQHEs)  without  Landau  levels  (58,  59). 

In  recent  years,  the  field  of  ultracold  atoms  has  made  re¬ 
markable  progress  in  engineering  topological  quantum  matter. 
An  artificial  gauge  field  (67)  has  been  realized  using  cold  atoms 
loaded  into  shaken  optical  lattices  (68,  69)  as  well  as  in  lattices 
with  laser-induced  tunneling  (55,  54,  70-72).  Various  topolog¬ 
ical  models,  including  Haldane’s  honeycomb  lattice  (73,  74), 
have  been  successfully  implemented.  Berry  curvature  and  to¬ 
pological  invariants  such  as  the  Chern  number  (74-80)  can  be 
measured.  Chiral  edge  currents  in  synthetic  quantum  Hall  lat¬ 
tices  are  also  observed  (81,  82).  Most  of  the  demonstrations 
so  far  focus  on  probing  topological  band  structures  and  single¬ 
particle  physics.  Realizing  strongly  interacting  topological 
phases  such  as  FQH  states,  however,  still  remains  elusive.  This 
in  part  is  due  to  limited  topological  bandwidth-to-gap  ratio,  but 
a  number  of  improved  schemes  (e.g.,  refs.  83  and  84)  have 
been  proposed. 

Coupling  cold  atoms  to  mobile  PCW  photons  also  allows  to¬ 
pological  band  engineering  and  band  flattening.  Moreover,  the 
pseudo  spin- 1/2  system  already  interacts  like  hard-core  bosons 
because  individual  atoms  that  participate  in  the  spin-exchange 
process  cannot  be  doubly  excited.  With  the  addition  of  tunable 
long-range  ZZ  interactions,  we  can  readily  build  many-body 
systems  that  should  exhibit,  for  example,  FQH  and  supersolid 
phases  (85),  providing  a  powerful  route  toward  realizing  strongly 
interacting  topological  phases. 

Chiral-Flux  Square  Lattice  Model.  The  first  example  discussed  here 
can  be  mapped  to  a  topological  flat-band  model  similarly  de¬ 
scribed  in  refs.  58,  59,  and  85.  The  topological  spin  Hamiltonian 
is  written  as  follows: 


#fiat  =Hq+H' 

HQ=tl^el^n±t2  ^2  aman  +  h.c 

(m,n)  ((m,n))  [15] 

H'  =  t3  ^2 

(((m,n))) 

in  which  we  define  0^=0^  and  am=o^s\  (.)  denotes  nearest 
neighbors  (NN),  and  t\  is  the  coupling  coefficient,  ((.))  [(((.)))] 
denotes  [next-]next-nearest  neighbors  (NNN)  [and  NNNN,  re¬ 
spectively]  with  t2  fe]  being  the  respective  coupling  coefficients. 
The  NN  coupling  phases  (pmn  =  ±</>  are  staggered  across  lattice 
sites,  where  the  phase  factor  </>  is  the  one  that  breaks  time  reversal 
symmetry  for  (p^0,ri7i  (with  n  G  Z).  Spin  exchange  between  next- 
nearest  neighbors  (NNN)  has  real  coefficients  ±t2  with  alternating 
sign  along  the  lattice  checkerboard  (Fig.  6).  One  can  show  that 


N 

Fig.  5.  Sideband  amplitude  for  Haldane-Shastry  model:  total  intensity 
(black)  J2\*a\2  ar|d  maximum/minimum  ratio  (red)  of  sideband  amplitudes 

\Xa\  as  a  function  of  N. 


Fig.  6.  Engineering  a  chiral-flux  square  lattice.  (/A)  Two  sublattices  ( nx  +  ny 
odd  or  even)  are  marked  by  blue  and  red  circles,  respectively.  Solid  lines 
mark  the  NN  hopping  with  phase  gain  cp  (arbitrarily  tuned)  along  the  di¬ 
rection  of  the  arrows.  Dashed  (dotted)  lines  mark  the  NNN  hopping  terms 
(coefficients  ±£2)-  NNNN  long-range  hopping  along  curved  lines  are  included 
to  assist  band  flattening.  Filled  arrows  indicate  the  propagation  of  pump 
electric  fields  sy  and  sXl  respectively;  see  text.  ( B )  Resulting  two-band  struc¬ 
ture  with  (t2,  t3,  <p)  =  (ti  /a/2,  ti  /4a/6,  n/4). 


already  Hq  has  a  small  bandwidth  with  nontrivial  Chern  number 
that,  choosing  t2=t\/y/2  and  (j)  =  n/ 4,  results  in  a  simple  band 
dispersion  Lo(k)=  ±s/2t\  a/3  +  cos (kx  4-  ky)cos(kx  - ky).  Adding 
H'  to  Hq  with,  for  example,  t3  =  1  /4v6*i  allows  us  to  engineer  an 
even  flatter  lower  band  whose  bandwidth  is  ~  1  %  of  the  band  gap. 

We  can  use  an  array  of  atoms  trapped  within  a  2D  PCW,  as 
in  Fig.  IB,  to  engineer  the  Hamiltonian  i/fiat  of  Eq.  15.  For 
simplicity,  we  assume  that  there  is  one  atom  per  site  although 
this  is  not  a  fundamental  assumption.8  As  shown  in  Fig.  6^4,  we 
need  to  engineer  spin  exchange  in  four  different  directions, 
namely,  x,y,x±y.  We  first  introduce  linear  Zeeman  shifts  by 
properly  choosing  a  magnetic  field  gradient  VB  (SI  Appendix  B: 
Proper  Choice  of  Ground- State  Energy  Shifts  in  2D  Models)  such 
that  8 a  =  \pBS7B  ■  Ara\,  where  ptB  is  the  magnetic  moment,  Ara 
are  vectors  associated  with  the  directions  of  spin  exchange: 


{Ara}„=xvaw,  =  {dx,dy,d{x+y),d(x-y)},  and  d  is  the  lattice 
constant.  To  activate  spin  exchange  along  these  directions  while 
suppressing  all  other  processes,  we  consider  a  simplest  case  by 
applying  a  strong  pump  field  of  amplitude  £2o  (frequency  coi)  to 
pair  with  sidebands  |£2a|  <C  |£2o|  of  detunings  cba  =  8a  to  satisfy  the 
resonant  conditions.  To  generate  the  desired  chiral-flux  lattice,  we 
need  to  carefully  consider  the  propagation  phases  ko  •  rn  (ka  •  rn) 
of  the  pump  field  (and  sidebands),  where  rn  =  d(nxx  +  nyy)  is  the 
site  coordinate  and  nxy  e  Z.  In  the  following,  we  pick  k  =  ka  =  n/d. 

We  can  generate  the  couplings  in  7/flat,  term  by  term,  as  follows. 
Staggered  NN  coupling  along  Ara=Xjy.  We  consider  the  strong  pump 
field  to  be  propagating  along y,  that  is,  Xo(rn)  =  (\Q.o\/2A)e~lny!l. 
At  the  NN  site  rm  =  rn  +  Arx,  it  can  pair  with  an  auxiliary 
sideband  of  detuning  a>x  =  8X  =  \pBV B  -  Arx\  with  Xx(rm)  = 
(|Qi  |/2 A)  [e~inyn  -  i^e~l^nx+1^]  to  generate  coupling  along  Ar*.  The 
sideband  is  formed  by  two  field  components  in  sy(t)  and  sx(t),  propa¬ 
gating  along  y  and  x,  respectively  (Fig.  6),  with  an  amplitude  ratio  of 
£  and  with  an  initial  n /2  phase  difference.  These  two  fields  are  used  to 
independently  control  real  and  imaginary  parts  of  the  spin-exchange 
coefficients.  Using  Eq.  8  under  the  condition  |£2o|  >  |£2«|,  the  cou¬ 
pling  rate  along  Axx  is  as  follows: 


Jm,n=JX0X*=t1 


yi+c1 


=  t\e 


,±i(p 


[16] 


where  t\  =J\Xo  \  \Xx  \  \J\  +  f2.  This  results  in  the  staggered  phase 
pattern  with  tunable  </>  =  tan-1  f .  The  NN  coupling  along  Aiy  can 


§ln  principle,  exact  physical  separations  between  trapped  atoms  do  not  play  a  significant 
role  with  photon-mediated  long-range  interactions.  One  may  also  engineer  the  spin 
Hamiltonian  based  on  atoms  sparsely  trapped  along  a  photonic  crystal,  even  without 
specific  ordering.  It  is  only  necessary  to  map  the  underlying  symmetry  and  dimensionality 
of  the  desired  spin  Hamiltonian  onto  the  physical  system. 
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Fig.  7.  Engineering  a  honeycomb-equivalent  topological  brick  wall  lattice. 
(A)  Unit  cell  of  a  honeycomb  lattice.  Solid  lines  mark  the  NN  hopping. 
Dashed  lines  mark  the  NNN  hopping  with  phase  gain  $  along  the  direction 
of  the  arrows.  ( B )  Unit  cell  of  a  brick  wall  lattice.  Solid  lines  indicate  the  NN 
hopping  as  in  A.  NNNN  hopping  (curved  dashed  lines)  and  NNN  hopping 
(diagonal  dashed  lines)  correspond  to  the  complex  NNN  hopping  in  A, 
making  the  two  models  topologically  equivalent.  (C)  Brick  wall  lattice.  Filled 
arrows  illustrate  the  pump  electric  fields.  (D)  Band  structure  of  the  brick  wall 
lattice,  plotted  with  cos^  =  3^3/43  (58). 


be  introduced  via  another  sideband  with  detuning  cby  =  dy  = 
\jtBVB- Ar_y|  andXy(rn  +  Ary)=  -(|£2i|/2A)  [e~l^y+^n +  il^e~inxn). 
NNN  coupling  along  Ara=xy,xy*.  The  sign  of  the  coefficient  depends 
on  the  sublattices.  To  engineer  these  couplings,  we  use  two 
sidebands  formed  by  field  components  in  sx(t),  with  detunings 
8a=sxyw‘  and  Xa(rn  +  Arw.)  =  ±(|Q2|/2A)e-,>r^+1)  at  NNN  sites. 
After  pairing  with  the  pump  field  Xq  at  site  rn,  the  resulting  ex¬ 
change  coefficients  are  Jm/l  =JX oX*=xyjy*  =  +t2(-l)nx~ny,  forming 
the  required  pattern  with  ^  =/|^b||^*y|- 

NNNN  coupling  along  2Ara=x,y.  We  use  two  sidebands  X^y  = 
|Q3|e_WI”y/2A,  propagating  along  y  with  detunings  2Sa=xy,,  to  in¬ 
troduce  the  real  coupling  coefficient  t$  =J\Xo\\X2x\. 

Summing  up,  all  of  the  components  in  the  Raman  field  can  be 
introduced  by  merely  two  pump  beams  propagating  along  x  and 
y  directions,  respectively.  In  SI  Appendix  C:  Pump  Field  Config¬ 
urations  for  Engineering  a  Chiral-Flux  Square  Lattice  Model ,  we 
explicitly  write  down  the  time-dependent  electric  field  that 
contains  all  of  the  sidebands. 

We  note  that  it  is  also  possible  to  simultaneously  introduce  both 
blue-detuned  (<5a>0)  and  red-detuned  ( S-a  =  -8a )  sidebands  in 
the  Raman  field  to  control  the  same  spin-exchange  term.  That  is, 
Jm,n  =J[Xo(rn)X*(rm)  +Xq  (rm)X_a(r„)],  which  has  contributions 
from  Xa  and  X-a  of  blue  and  red  sidebands,  respectively.  Arranging 
both  sidebands  with  equal  amplitudes  lead  to  equal  contribu¬ 
tions  in  the  engineered  coupling  coefficient.  This  corresponds 
to  applying  amplitude  modulations  in  the  pump  electric  field. 
In  real  experiments,  amplitude  modulation  can  be  achieved 
by,  for  example,  the  combination  of  acoustic-optical  modulators, 
and  optical  IQ-modulators. 


NNNN  (along  2Ary)  couplings  with  the  same  strength  and 
with  a  coupling  phase  (frmn  =  ±</>,  which  alternates  sign  across 
two  sublattices.  Thus,  our  target  Hamiltonian  is  given  by  the 
following: 


H  —  t\  (plnan  +h.C.)  +t2  {el^mn +h.c.),  [17] 

(m,n)  {m,n} 


where  (. )  denotes  NN  pairs  in  the  brick  wall  configuration  (Fig. 
7)  and  t\  is  the  coupling  coefficient.  Note  that,  for  simplicity,  we 
discuss  a  special  case  where  all  NN-coupling  coefficients  from  a 
brick  wall  vertex  are  identical.  The  second  summation  in  Eq.  17 
runs  over  both  NNN  and  NNNN  pairs  with  identical  coupling 
coefficient  ^  and  alternating  phase  (j)mn  =  ±(p  (Fig.  7). 

As  in  the  previous  case,  we  use  a  strong  pump  field  (propa¬ 
gating  along  y),  as  well  as  several  other  weak  sidebands  to  gen¬ 
erate  all  necessary  spin-exchange  terms.  Detailed  descriptions  on 
engineering  individual  terms  can  be  found  in  SI  Appendix  D: 
Pump  Field  Configurations  for  Engineering  a  Topological  Spin 
Model  in  a  Brick  Wall  Lattice.  The  most  important  ingredient, 
discussed  here,  is  that  we  can  generate  checkerboard-like  NN 
coupling  (along  x),  with  /m>„  =JXqX*  =  {t\  /2)  [1  -  (-l)"*-"*'].  This  is 
achieved  by  using  a  sideband  of  detuning  8X  and  amplitude 
Xx  =  (|Q|/4A)  \e~inyn  +  ^e~l^nx+V)n\  at  position  rm  =  rn  +  Arx,  formed 
by  two  fields  propagating  along  y  and  x,  respectively.  If  both 
fields  have  the  same  amplitude  (f  =  1),  they  either  add  up  or 
cancel  completely  depending  on  whether  nx  -  ny  is  odd  or  even. 
If  one  applies  the  same  trick  toward  NN  coupling  along  y,  but 
with  f^l,  the  coupling  amplitude  modulates  spatially  in  a 
checkerboard  pattern.  Essentially,  all  three  NN  terms  around  a 
brick  wall  vertex  can  be  independently  controlled,  opening  up 
further  possibilities  to  engineer,  for  example,  Kitaev’s  honey¬ 
comb  lattice  model  (87,  88). 

For  physical  implementations,  again  only  two  pump  beams  can 
introduce  all  components  required  in  the  Raman  field,  which  is 
very  similar  to  the  previous  case.  We  stress  that,  by  merely 
changing  the  way  the  Raman  field  is  modulated,  one  can  dy¬ 
namically  adjust  the  engineered  spin  Hamiltonians  and  even 
the  topology,  as  we  compare  both  cases.  This  is  a  unique  feature 
enabled  by  our  capability  to  fully  engineer  long-range  spin 
interactions. 

Moreover,  many  of  the  tricks  discussed  above  can  also  be 
implemented  in  ID  PCWs.  It  is  even  possible  to  engineer  a 
topological  ID  spin  chain,  by  exploiting  long-range  interac¬ 
tions  to  map  out  nontrivial  connection  between  spins.  For 
example,  our  method  can  readily  serve  as  an  realistic  ap¬ 
proach  to  realize  a  topological  ID  spin  chain  as  recently 
proposed  in  ref.  89. 

XXZ  Spin  Hamiltonian  with  Tunable  Interaction  1  /r'K  In  the  last  ex¬ 
ample,  we  highlight  the  possibility  of  engineering  a  large  class  of 
XXZ  spin  Hamiltonians,  which  were  studied  extensively  in  the 
literature  because  of  the  emergence  of  frustration  related  phe¬ 
nomena  (60,  90-96)  and  their  intriguing  nonequilibrium  dynamics 
(45-50).  An  XXZ  Hamiltonian  is  typically  written  as  follows: 

Hxxz  =  -B  <  +  E  fr~  [co s(0K<4  +  sin(0)  («,  +  «,)}  > 

n  n<m  rn’m 

[18] 


"Honeycomb"-Equivalent  Topological  Lattice  Model.  To  further 
demonstrate  the  flexibility  of  the  proposed  platform,  we  create 
Haldane’s  honeycomb  model  (73)  via  a  topologically  equivalent 
brick  wall  lattice  (74,  86).  Here,  we  engineer  the  brick  wall 
configuration  using  the  identical  atom-PCW  platform  discussed 
in  the  previous  example.  Mapping  between  the  two  models  is 
illustrated  in  Fig.  7  A  and  B ,  which  contains  the  following  two 
nontrivial  steps:  (/)  generating  a  checkerboard-like  NN-exchange 
pattern  in  the  x  direction;  ( ii )  obtaining  NNN  (along  Ar^,*)  and 


where  an  effective  magnetic  field  B  controls  the  number  of  excita¬ 
tions,  rm,m  =  |rn  -  rm|,  and  the  parameter  0  determines  the  relative 
strength  between  the  ZZ  and  XY  interactions.  This  class  of  spin 
models  has  been  previously  studied,  but  mostly  restricted  to  nearest 
neighbors  (90-93)  or  dipolar  (rj  =  3)  interactions  (60,  94,  95). 

In  our  setup,  one  can  simulate  XXZ  models  with  arbitrary  rj  by 
first  introducing  unique  ground-state  energy  shifts  at  each  of  the 
separation  rn  -  rm,  and  then  applying  a  strong  pump  field  of  am¬ 
plitude  together  with  Nd  auxiliary  fields  Qa  of  different  detunings 
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XX Z  with  4x4  spin  lattice 
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Fig.  8.  Mean  magnetization  M/N  for  a  system  with  N=  16atoms  in  a  square 
lattice,  restricted  to  A/exc<8  excitations  and  rj=  1  (A),  2  ( B ),  3  (C),  and  NN 
couplings  (D). 


of  the  many-body  states  may  be  probed  via  guided  photons  in  an¬ 
other  propagation  mode  along  the  PCW  (98). 

Limitations  and  Error  Analysis.  Until  now,  we  have  mainly  focused 
on  how  to  engineer  Hq  in  an  ideal  situation.  We  neglected  spon¬ 
taneous  emission  or  GM  photon  losses  and  considered  that  the 
energy  gradient  (or  8,  the  ground-state  energy  difference  between 
nearest  neighboring  atoms)  can  be  made  very  large  compared 
with  the  interaction  energy  scales  that  we  want  to  simulate 
(|<5|  >  \Jm,m+i !)•  Because  the  effect  of  finite  cooperativities  was 
considered  in  detail  in  refs.  36  and  37,  and  their  conclusions 
translate  immediately  to  our  extension  to  multifrequency  pumps, 
in  this  work  we  mainly  focus  on  the  effect  of  finite  8.  In  addition, 
we  also  discuss  the  effects  of  AC  Stark  shifts  as  in  Eq.  9,  and  its 
error  contributions,  together  with  other  possible  error  sources. 

Corrections  Introduced  from  Higher  Harmonics:  A  Floquet  Analysis. 

We  discuss  errors  and  the  associated  error  reduction  scheme 
following  a  Floquet  analysis  with  multifrequency  driving  (99, 
100),  applicable  mainly  to  ID  models.  Including  all  of  the  time- 
dependent  terms  in  a  multifrequency  pumping  scheme,  we  have 
(Eq.  5)  H(t)  =  Y  HpeipSt,  where  Hp  represents  the  part  that  oscil¬ 
lates  at  frequency  pS.  This  Hamiltonian  has  a  period  T  =  2 n/8.  It 
can  be  shown  that  at  integer  multiples  of  T,  the  observed  system 
should  behave  as  if  it  is  evolving  under  an  effective  Hamiltonian#: 


to  introduce  spin  interactions  at  each  separations.11  Moreover,  the 
parameter  6  that  determines  the  ratio  between  ZZ  and  XY  in¬ 
teraction  can  be  controlled  by  using  different  pump  intensities  in  the 
stroboscopic  steps  (SI  Appendix  E:  PCW  and  Pump  Field  Configu¬ 
rations  for  Engineering  an  XXZ  Spin  Hamiltonian  with  1/ rn 
Interaction). 

To  illustrate  physics  that  can  emerge  in  the  first  experimental 
setups  with  only  a  few  atoms,  we  study  the  total  magnetization  of  a 
small  square  lattice  of  ns  x  ns  (=N)  16  atomic  spins.  We  apply  exact 
diagonalization  restricting  to  Nexc  <  8  excitations  for  N  =  16  spins 
and  cover  one-half  of  the  phase  diagram  with  B  >  0.  In  Fig.  8,  we 
explore  the  mean  magnetization  of  the  system  M/N  =  \  Y?  (of)  /N  as 
a  function  of  B  and  6  for  rj  =  1  (A),  2  ( B ),  3  (C),  and  NN  cou¬ 
plings  (D).  At  6  =  0,  the  system  behaves  classically  showing  the  so- 
called  “devil  staircase”  (97)  of  insulating  states  with  different 
rational  filling  factors  and  crystalline  structures.  As  already  explored 
in  ref.  95  for  ID  dipolar  systems  (rj  =  3),  the  presence  of  long-range 
interactions,  compared  with  nearest-neighbor  models,  lead  to 
stronger  frustration  effects.  This  manifests  in  the  magnetization 
diagram  with  an  asymmetry  between  6  $  0.  In  Fig.  8,  we  show  that 
longer-range  interactions  lead  to  even  higher  degree  of  asymmetry. 
Moreover,  in  refs.  60  and  95,  it  was  discussed  that  long-range 
coupling  leads  to  the  formation  of  supersolid  phases,  in  which 
crystalline  structure  and  long-range  order  coexist.  These  may  be 
even  more  favored  by  longer-range  interactions.  Full  character¬ 
ization  of  the  phases  diagram  is,  however,  beyond  the  scope  of  this 
paper  and  will  be  discussed  elsewhere. 

Another  especially  interesting  arena  is  the  behavior  of  strongly 
long-range  interacting  systems  (rj  smaller  than  the  lattice  di¬ 
mension  D )  under  nonequilibrium  dynamics.  It  has  recently  been 
predicted  to  yield  “instantaneous”  transmission  of  correlations 
after  a  local  quench  (45,  47,  48,  96),  breaking  the  so-called  Fieb- 
Robinson  bound. 

Finally,  it  is  interesting  to  point  out  that  magnetization  can  be 
measured  by  first  freezing  the  interaction  (via  shutting  off  the  pump 
lasers)  followed  by  atom  number  counting  using  state-dependent 
fluorescence  imaging.  Coherence  and  off-diagonal  long-range  orders 


#eff,l 


[HP,H-p] 

P 


1  v[[j Hp,Ho],H-p\  +  [[H-p,Ho\,Hp ] 
2  S2^  P2 


[19] 


This  means  that  the  leading  error  in  our  simple  scheme  would  be 
on  the  order  of  J2/8,  where  J  is  the  simulated  interaction 
strength.  However,  we  note  that  if  Hp  =  ±H-P,  the  leading  error 

term  Y[Hp>H-p\/ (P&)  should  vanish.  In  other  words,  first-order 

error  vanishes  if  Hp  is  either  symmetric  or  antisymmetric  under  a 
time  reversal  operation  X.  Although  the  original  Hamiltonian 
H(t)  does  not  necessarily  possess  such  symmetry,  it  is  possible 
to  introduce  a  two-step  periodic  operation  H2^&v  =  {H,JrH, 
H,SEH,  . . . }  to  cancel  the  first-order  error  while  keeping  the 
time-independent  part  //2step,o=#o  identical.  This  results  in  an 
effective  Hamiltonian  in  the  Floquet  picture: 

He ff,2  =H0  +  Hqu2  «H0+I  Y.  Hf  i&lEAA,  [20] 

8  z p 
p  1 

where  Hp  is  the  (operator)  Fourier  coefficient  of  the  two-step 
Hamiltonian  and  the  leading  error  reduces  to  the  order  of  J3  /82. 

To  achieve  the  time  reversal  operation,  we  must  reverse  the 
phase  of  the  driving  lasers,  as  well  as  the  sign  of  the  energy 
offsets  between  the  atoms.  Specifically,  we  can  engineer  a  peri¬ 
odic  two-step  Hamiltonian  by  first  making  the  system  evolve 
under  presumed  Ho  (along  with  other  time-dependent  terms)  for 
a  time  interval  T,  and  then,  for  the  next  time  interval  T,  we  flip 
the  sign  of  the  energy  gradient,  followed  by  reversing  the  prop¬ 
agation  direction  of  the  Raman  fields  such  thatXa  -+X*  in  Eq.  5. 
As  a  result,  all  of  the  time-dependent  Hamiltonians  Hp ,  \/p  ±  0, 
become  H_p  in  the  second  step,  resulting  in  Hp  =  (-1  ffH~p  re¬ 
quired  for  error  reduction;  whereas  the  time-independent 
Hamiltonian  Ho  remains  identical  in  the  two-step  Hamiltonian. 
See  SI  Appendix  F:  Error  Reduction  and  Analysis  for  more  discussions. 


^To  simulate  a  square  lattice  of  nsxns  ( =N)  atomic  spins,  we  find  that  the  number  of 
different  distances  grows  as  A/c/  =  (ns (ns  + 1 )  —  2) /2#  which  is  linearly  proportional  to  the 
number  of  atoms  Ndoc.N. 


#When  the  measurement  time  is  incommensurate  with  period  T,  small-amplitude  and 
fast-oscillating  spin-dynamics  due  to  time-dependent  terms  in  Eq.  5  manifest  as  extra 
errors;  see  the  discussion  about  micromotion  in  refs.  99  and  100. 
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Fig.  9.  (/A)  Comparison  of  ground-state  energy  error  |(£0 -£/)/£0|  and  ( B ) 

ground-state  overlap  |('P0|'P/)|  as  a  function  of  N  for  Hamiltonians  Heff(1 
(black)  and  He ff(2  (red)  with  detuning  bjp  =40. 


Numerical  Analysis  on  the  Haldane-Shastry  Spin  Chain.  We  now 

analyze  numerically  and  discuss  error  on  one  particular  example. 
For  numerical  simplicity,  we  choose  the  Haldane-Shastry  model 
as  its  ID  character  makes  it  numerically  more  accessible.  How¬ 
ever,  the  conclusions  regarding  the  estimation  of  errors  can  be 
mostly  extended  to  other  models.  As  we  have  shown  in  Eq.  13, 
the  Haldane-Shastry  Hamiltonian  is  composed  by  an  XY  term 
plus  a  ZZ  term  that  we  can  simulate  stroboscopically.  As  we 
already  analyzed  the  Trotter  error  due  to  the  stroboscopic  evo¬ 
lution,  here  we  focus  on  the  XY  part  of  the  Hamiltonian,  which 
reads  as  follows: 


#«*=££ -t 


Jo 


Si  sin  (nn/N) 


[21] 


where  Jo  =Jn2/N2.  Following  the  prescribed  engineering  steps, 
the  total  time-dependent  Hamiltonian  resulting  from  multiple 
sidebands  can  be  written  as  H(t)  =  Y Hpeip5t ,  with  the  following: 

p 

N  N-m 

hp = £  £  (4(P)«+n  l22] 

m— 1  n— 1 

and  we  have  defined  J„,{p]  =  Y,aJ)LXxp8n-p,/i-a-  Here,  Xa  are 
fixed  such  that  Ho  =Hus,xy 

To  illustrate  the  effect  of  error  cancellations,  we  consider  first 
a  scenario  where  we  directly  apply  Eq.  22.  To  leading  order,  the 
effective  Hamiltonian  is  as  in  Eq.  19.  We  then  analyze  the 
two-step  driving,  using  the  effective  Hamiltonian  //eff,2  in  Eq.  20, 
with  Hp  given  by  Eqs.  S45  and  S46. 

We  calculate  the  ground-state  energies  and  eigenvectors  of 
Ho ,  He ff;i,  and  He ff;2,  which  we  denote  as  E0)  1)2  and  |T0jij2),  re¬ 
spectively,  for  different  number  of  atoms  and  different  ratios  of 
S/Jo.  The  results  are  shown  in  Figs.  9  and  10.  In  panels  A ,  we 
show  the  error  in  absolute  value  with  respect  to  the  ideal 
Hamiltonian  Ho.  Interestingly,  due  to  particular  structure  of  |To) 
and  Hp ,  one  can  show  that  (T0|[F/p,//_;7]|To)  and  the  first- 
order  correction  to  the  energy  vanishes.  This  is  confirmed  in  Fig. 
10 L4,  where  we  found  that  the  error  actually  scales  with  1  /S2. 

Moreover,  it  is  also  enlightening  to  compare  the  overlap  of  the 
ground  states  as  shown  in  Figs.  9 B  and  10 B.  We  only  compute 


the  even-atom  number  configuration  as  the  odd  ones  are  de¬ 
generate  and  therefore  the  ground  state  is  not  uniquely  defined. 
We  see  that  the  ground-state  overlap  of  //eff,2  is  several  orders 
of  magnitude  better  than  the  one  with  Moreover,  its 

dependence  on  S  is  better  than  the  1/S2  expectation. 

The  Role  of  Time-Dependent  Stark  Shifts  in  the  Error  Analysis.  In  the 

previous  discussions,  we  have  dropped  the  contribution  of  the 
time-dependent  Stark  shifts: 


mp-l 

tfae(0  =  -££  * 

n  a>j3 


2A 


[23] 


where  coajj  =  cba-  In  SI  Appendix  F:  Error  Reduction  and 
Analysis ,  we  discuss  its  role  in  the  effective  Hamiltonian,  using 
the  Floquet  error  analysis.  To  summarize,  we  evaluated  the  error 
in  the  two-step  driving  scheme  in  various  configurations. 

Generic  Hamiltonians  with  translational  invariance.  By  translational 
invariance,  we  mean  that  there  are  no  site-dependent  spin  in¬ 
teractions,  and  the  spin-exchange  coefficients  remain  identical  as 
we  offset  the  spin  index  by  one  or  more.  This  means  that  all 
components  in  the  pump  field  should  drive  the  system  with 
uniform  optical  phases  as  in  the  Haldane-Shastry  model  dis¬ 
cussed  above.  The  error  by  Hac{t)  averages  out  to  zero  in  the 
Floquet  picture.  In  the  butterfly  scheme,  however,  both  \g)  and 
| s)  states  are  pumped  and  they  may  be  shifted  differently.  This 
leads  to  slight  modifications  in  the  engineered  XX  and  YY  terms. 
Models  containing  sublattices.  For  topological  models  that  contain 
sublattices,  as  in  our  examples,  the  pump  fields  are  not  perfectly 
transverse  and  Stark  shifts  are  site  dependent,  resulting  in  non¬ 
vanishing  error.  For  realistic  PCW  realizations,  one  should  set 
moderate  pump  detuning  7/ A  >  0(1)  such  that  leading  error 
contribution  will  be  <J3/S2,  and  for  <5»/,  the  Stark  shift 
terms  may  be  ignored. 

Stark  shift-dominated  regime.  It  may  be  possible  that  our  sublattice 
models  be  purposely  driven  with  large-amplitude  pumps  such 
that  |£2|2/A  >  S.  Stark  shift  contributions  would  become  important 
in  the  resulting  spin  dynamics.  However,  if  we  choose  a  large  pump 


S/Ji 


Fig.  10.  (A)  Comparison  of  ground-state  energy  error  |(£0-£,)/£o|  and  ( B ) 
ground-state  overlap  |('P0|'P/)|  as  a  function  of  5/J ^  for  Hamiltonians  Heff,i 
(black)  and  Hef f  2  (red)  for  A/  =  12  atoms. 
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detuning  A  >7,  the  dominant  error  contribution  can  in  fact  be 
written  in  the  following  simple  form: 

Herra  *  J21 2 3 4 5 6 7 8 9 10 11 12m,n  +  **■)  -  [24] 

m,n 

where  Am,n  is  a  site-dependent  amplitude.  In  a  special  case  that 
only  two  sublattices  are  present,  as  in  our  examples,  we  note  that 
Am,n  may  only  depend  on  the  distance  and  is  site  independent. 
This  “error”  term  would  then  uniformly  modify  the  XY  coupling 
strengths  to  a  new  value: 

J m,n  =  (l  m,n)Jm,ri'  [25] 

The  next  leading  order  errors  are  a  factor  of  ~7/A  smaller  than 
this  leading  Stark  shift  contribution,  suggesting  we  can  always 
increase  the  detuning  A,  while  keeping  |£2|/A  constant,  to  reduce 
the  error  contribution. 

Other  Error  Sources  and  Heating  Effects.  Apart  from  errors  arising 
from  multifrequency  driving,  there  are  other  common  error 
sources  in  cold  atoms  that  we  have  not  considered  so  far,  such  as 
motional  heating.  In  the  PCW  platform,  atoms  are  tightly  con¬ 
fined  with  a  trap  depth  more  than  three  orders  of  magnitude 
larger  than  the  recoil  energy,  rendering  well-separated  motional 
bands  such  that  effects  like  interband  heating  (101)  can  be 
suppressed.  Spin-exchange  rates  in  the  PCW  platform,  however, 
can  be  adjusted  to  1  MHz  >  \Jxy*\  >  1  kHz  so  that  the  many- 
body  time  scales  (<^C  1  ms)  can  be  much  faster  than  those  asso¬ 
ciated  with  motional  heating. 

In  fact,  spin  temperature  can  be  decoupled  from  real  atomic 
temperature  while  simulating  the  spin  models.  For  example,  one 
can  polarize  atomic  spins  initially  in  a  strong  magnetic  field 
( B  >  to  approximate  a  zero-temperature  paramagnetic 

phase  (17).  The  magnetic  field  can  then  be  ramped  down  adia- 
batically  to  the  final  value  of  the  desired  spin  model.  Limitations 
to  adiabaticity  and,  therefore,  to  the  accessible  spin  temperature 
will  ultimately  be  limited  by  the  fidelity  of  the  spin-exchange  (36, 
37)  or  by  motional  heating  that  leads  to  dephasing,  whichever 
gives  a  more  stringent  bound. 

Conclusions  and  Outlook 

In  this  paper,  we  have  shown  that  atom-nanophotonic  systems 
present  appealing  platforms  to  engineer  many-body  quantum 
matter  by  using  low- dimensional  photons  to  mediate  interaction 
between  distant  atom  pairs.  We  have  shown  that,  by  introducing 
energy  gradients  in  ID  and  2D,  and  by  applying  multifrequency 


Raman  addressing  beams,  it  is  possible  to  engineer  a  large  class 
of  many-body  Hamiltonians.  In  particular,  by  carefully  arranging 
the  propagation  phases  of  Raman  beams,  it  is  possible  to  introduce 
geometric  phases  into  the  spin  system,  thereby  realizing  nontrivial 
topological  models  with  long-range  spin-spin  interactions. 

Another  appealing  feature  of  our  platform  is  the  possibility  of 
engineering  periodic  boundary  conditions,  as  explicitly  shown  in 
the  ID  Haldane-Shastry  model,  or  other  global  lattice  topology  by 
introducing  long-range  interactions  between  spins  located  at  the 
boundaries  of  a  finite  system.  Using  2D  PCWs,  for  example,  it  is 
possible  to  create  previously  unavailable  spin-lattice  geometries  such 
as  Mobius  strip,  torus,  or  lattice  models  with  singular  curvatures  such 
as  conic  geometries  (102)  that  may  lead  to  localized  topological 
states  with  potential  applications  in  quantum  computations. 

We  emphasize  that  all  of  the  pairwise- tunable  interactions  can 
be  dynamically  tuned  via,  for  example,  electro-optical  modula¬ 
tors  at  timescales  much  faster  than  that  of  characteristic  spin 
interactions.  Therefore,  the  spin  interactions  can  either  be  adi- 
abatically  adjusted  to  transform  between  spin  models  or  even  be 
suddenly  quenched  down  to  zero  by  removing  all  or  part  of  the 
Raman  coupling  beams.  We  may  monitor  spin  dynamics  with 
great  detail:  after  we  initially  prepare  the  atomic  spins  in  a 
known  state  by,  say,  individual  or  collective  microwave  address¬ 
ing,  we  can  set  the  system  to  evolve  under  a  designated  spin 
Hamiltonian,  followed  by  removing  all  of  the  interactions  to 
“freeze”  the  dynamics  for  atomic  state  detection.  Potentially,  this 
allows  for  detailed  studies  on  quantum  dynamics  of  long-range, 
strongly  interacting  spin  systems  that  are  driven  out-of-equilibrium. 
The  dynamics  may  be  even  richer  because  the  spins  are  weakly 
coupled  to  a  structured  environment  via  photon  dissipations.  We 
expect  such  a  platform  may  bring  novel  opportunities  to  the  study  of 
quantum  thermalization  in  long-range  many-body  systems,  or  for 
further  understanding  of  information  propagation  in  a  long-range 
quantum  network. 
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